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I.  Formulation  of  the  problem.. 

In  geophysical  prospecting  and  seismology  there  are  many  problems  which  are 
governed  by  the  equation: 

Au(x)  +  uj^n{z)u{x)  =  -*(a;)  (1) 

where  x  =  (x^is)  G  R'  ,  =  (ii.xj)  €R^  z  =  xs  <  0 ,  uj  >  0 .  In  acoustic 
prospecting,  S{x)  is  a  point  source  situated  at  the  origin,  u(i)  is  the  pressure, 
w  is  the  wave  frequency,  and  c{z)  =  n"^/'(z)  is  the  velocity  of  the  wave  at  the 
depth  z  .  The  inverse  problem  is  to  find  n(z)  from  the  given  data:  u(x,a;)  ,  where 
X  E  P  :=  {x  :  X3  =  0}  and  0  <  a;  <  wq  ,  ojq  is  a  fixed  small  frequency.  Let  us 
give  a  sketch  of  our  approjich  following  [l]. 
Fourier  transform  (1)  in  x^  to  get 

u"-X^u  +  q{z)u=-S{z), 

where 


u' 

u(A,z) 
A 

Assume  that  A'  >  go  :=      max,  \n{z)\ ,  then  the  integral  equation 

ti(A,  z)  =  g{X,  z)  +  g{X,  z  -  z')q{z')u{X,  z')dz'  (2) 

J  — OO 

where  g{X,z)  =  ,  is  uniquely  solvable  by  iterations,  and  u  is  analytic  in  A 

in  the  region  ReX  >  q^^'  .  Note  that 

u{X,z)  =  g{X,z) +oj^  f      giX,z-z')n{z')g{z')dz'  +  0{u*),as    u  ^  0. 

J  —  OO 

Therefore 

/(z.A)=  lim(u-,)a;-^  =  y_^  '--—n{zY-^dz'. 
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=  du/dz, 

/ +  00 
e-^^  '  u{x\z)dx^, 
-OO 

=  (Ai,A2), 

=  xl  +  xl 


This  low  frequency  limit  was  used  in  [2],  [3].  Set  z  =  0 ,  then 

Fix)  :=  4A V(0,  A)  =  e-'^\'ln{z)dz,    A  >  0.  (3) 

Since  the  data  are  given  on  the  plane  P  :=  {x  :  xs  =  0}  ,  u(A,0)  and  F{X) 
are  known.  Therefore  the  inverse  problem  is  reduced  to  finding  n{z)  from  the 
knowledge  of  F{X) . 

Assume  that  n(z)  =  1  for  z  >  0  and  n(z)  =  no  =  const  for  z  <  -rf,  where 
d  is  a  certjiin  depth.  Then  (3)  becomes 

l\-''\M^)dz  =  rPiX)  (4) 

where  il>{X)  =  F{X)  -  ^  -  .  A  change  of  variables,  2A  =  p ,  z  =  -i , 

transforms  (4)  into 

/  e-P*h{t)dt  =  ^p),  p>0 
Jo 

where  h{t)  =  n{-t) ,  =  V(f ) .   The  function  h{t)  can  be  found  by  the 

method  given  in  [3],  [4].  The  problem  of  inverting  the  Laplace  transform  of  a 
compactly  supported  function  from  the  real  axis  is  solved  analytically  in  [5]. 

We  assume  that  n(z)  =  1  for  z  >  0  (above  the  ground),  n(z)  =  no  ,  for 
z  <  -d,  no  >  0,  d>  0,  and  n(z)  =  ny  for  zy  <  z  <  zy+i  ,  1  <  j  <  m ,  zi  = 
^m+i  =  0 .  So  there  are  m  homogeneous  layers  in  the  region  —    <  z  <  0  . 
Therefore  (4)  yields 

m 

^n,[e-'^\''^^\-e-''M]  =  <f>{X):=2XfP{X),    A  >  0.  (5) 

3  =  1 

This  problem  is  ill-posed,  it  is  very  sensitive  to  noise  in  the  data.  One  wishes 
to  find  an  efficient  way  to  recover  parameters  {ny}^i  ,  {z,}^^i  from  the  knowl- 
edge of  ^(A) .  In  this  paper,  we  give  a  numerical  solution  to  this  problem  and 
demonstrate  its  efficiency.  We  assume  in  the  beginning  that  the  number  m  of 
layers  is  known.  In  section  IV  a  method  is  given  for  estimating  this  number  from 
the  data. 
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n.  The  description  of  the  algorithm. 

Let  {<^{Xi)}iLi  be  given,  where  {Ai}^i  is  a  set  of  positive  numbers.  Re- 
placing A  by  Xi  for  »  =  1, 2, . . . ,  iV  in  (5)  yields 

m 

.[e-2A,>,+.|  _  ^-2A,|,,||  ^  ^(^.j^      =  1,2,. . .  (6) 

3=1 

Tliis  is  a  nonlinear  system  for  {n,}^!  and  {zj}^^i  .  It  can  be  solved  by  an 
optimization  method.  Let  us  describe  the  method.  Set 

M^u  tfa,  •  •  • ,  Wm,  yi,  yj, . .  • ,  ym)  :=  J2    [«~'^'*'+'  -  e"'-^**^], 

3=1 
N 

G(«;i,u;,,...,«;^,yi,yj,...,y^)  :=       /?(tx;i,  ti;,,  yj,  yj, y„), 

j=i 

/<  =  ^-?i(A0,«  =  l,2,...,iV. 

We  want  to  find  {ny}^i  and  {zj}^^^  by  minimi2ing  G ,  the  target  function. 
Note  that  here  the  variables  wi,  u;,, . . . ,  u;„,  yi,  yj, . . . ,  y„  are  not  free  variables. 
They  satisfy  the  constraints: 

gj  :=  Wj  >  0, 1 
Ji+m  :=  yy  >  0,  >  J  =  l,2,...,m.  (7) 
93+2^  :=  yy+i  -  y,-  >  0,  J 

Introduce  parameters  Mj  ,  j  =  1, 2, . . . ,  3m  ,  then  the  problem  is  reduced  to  the 
following  non-constrained  optimization  problem: 

■f(«'i,wj,...,u;„,yi,y2,...,y^)  =  min 

where  F  =  G  +  ESf+i // ,  fi+s  =  M,|  min(0,50P  ,  i  =  1,2,. ..,3m,  are 
exterior  penalty  functions  [6],  and  Mi,  i  =  1,2, ...  ,3m ,  are  given  large  positive 
numbers  chosen  so  that  A  =  0 ,  t  =  1, 2, . . . ,  3m  ,  at  the  minimizer  of  F  ;  therefore 
conditions  (7)  are  satisfied  at  this  minimizer. 
Newton's  method  for  minimizing  F{x)  is 

X=  (xi,X2,...,Z2Ta)  =  (wi,ttf2,...,tWm,yi,y2,.-.,yn»). 
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Here  VF{x)  is  the  gradient  of     ,  and  F"  is  the  matrix  ( ^^^)  ,  i  =  1,2 
2"!  ,  J  =  1,2, . . .  ,2m  .  Newton's  method  has  qujidratic  convergence  necir  a  local 
minimizer.  But  it  requires  calculation  of  second  derivatives  of  the  target  function. 
One  can  estimate  the  second  derivatives  of  F{x)  by  the  first  derivatives  of  /j(x) , 
i  =  1,2,... ,M  :=  N  +  Zm. 

Let  be  an  approximation  of  the  minimizer  x*  of  F .    To  construct 

x(A:+i)  ^  we  approximate  /^(x)  by  the  Unear  functions 

3=1 

•^Vi(a:)=  -^,.  =  l,2,...,M,y=l,2,...,2m. 

So  in  place  of  F{x) ,  one  minimizes  Ei^i(^,-*^(«))'  •  Let  us  denote  this  function 
by  the  same  letter  F[x)  .  Then 

=  2(J(*))''(J(*)), 

=  ( Jo (xW)),  t  =  1, 2, . . . ,  M,y  =  1, 2, . . . ,  2m 
Vii'(x(*))  =  2(j(*))''/(i(*)) 

/(x)= 

V/m(x) 

Thus  the  iterative  procedure  is 

^(fc+i)  =       _  ((j(*))r(j(fc)))-i(j(*))T^(^(fc)) 

This  scheme  is  easy  to  implement  and  our  results  show  its  eflRciency.  Note  that 
is  always  ill-conditioned  when  i^*)  is  sufficiently  close  to  a  mini- 
mizer, and  sometimes  is  ill-conditioned  even  when  x^*)  is  not  close 
to  a  minimizer.  To  overcome  this,  we  change  the  search  direction  a  little,  so  the 
modified  procedure  is 

x(*+»)  =  x(*)  -  [(JW)''(J(*))  +  a,/]-Hj(*))''/(x(*))  (9) 

where  ajt  >  0  is  a  parameter.  The  choice  of  a*  is  described  in  the  outline  of  the 
program  below;  is  stored  in  the  cell  named  a .  This  is  Marquardt's  method 
which  is  a  combination  of  the  least  squares  method  and  the  gradient  method. 
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Let  us  outline  the  zdgorithm: 

One  starts  with  an  initial  point  a:(°)  ,  the  data  {Ai}^i  ,  {^(A,)},^i ,  m  ,  a 
positive  integer,  e ,  a  smeJl  positive  number,  d ,  a  positive  number,  parameters 
ao  >  0 ,  u  >  0,  {Afy}y^i  ,  each  Mj  is  a  positive  large  number.  IW  is  a  count 
number,  the  meaning  of  which  can  be  seen  clearly  from  the  following  outline. 

Step  1  :  0=^  k,  oq  =>  a  . 

Step  2  :  -1  =>  iW  ,  a/v  =>•  a. ,  calculate: 

(/i(x(*)),/,(xW)....,/M(xW)r, 

M 

»=i 

/  a/i(x'*')    3/i(x'*')      a/i(x'*')  \ 

9X1  9X2  dX3n> 

a/3(x'*')    a/2(x'*')  a/2(x"')) 

9xi  3x2  dX2m 


a/A^(x<*')    a/M(xt*))  a/M(xt*)) 

9X1  9X2  9x2m 

x(*'+i)  =  xW  +  P,    and  F{x^% 

where  aj^  is  a  real  number  in  cell  a . 

Step  4  :  If  is  less  than  F{x^''^) ,  then  pass  to  the  sixth  step, 

otherwise  do  the  next  step. 

Step  5:  IW  +  l=i^  IW  ,  au=>a.  U  \\  (j(*))''/(a;(*))  ||<  e ,  then  pass  to 
the  seventh  step,  otherwise  go  hack  to  the  third  step. 

Step  5  :  If  II  (J(*))^/(x(*))  ||<  e ,  then  go  to  the  seventh  step;  otherwise,  let 
k  +        k,  =►  x(*)  ,  and  go  back  to  the  second  step. 

Step  7  :  The  execution  stops. 

Since  there  is  no  guarantee  that  the  minimizer  found  by  this  algorithm  is  a 
global  one,  we  demand  that  the  optimal  point  x*  should  satisfy  two  inequalities 
II  ||<  c  and  F{x*)  <  ei  to  be  a  global  minimizer.  Here  ei  >  0  is 

a  small  number  and  c  >  0  is  the  number  given  at  the  start. 


/(x(*))  = 
f  (x(*))  = 

J(*)  = 

step  S  :  Calculate: 
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ni.  Nmnerical  results  and  practical  reconunendations. 
Rewrite  (5)  as 

n„»  +  (n^_i  -  n^)e-2^'«  +  (n^.j  -  rz„»_i)e-'^'~-'  +  •  •  • 
+  (ni  -  n2)e-^^'^  -  m*"'^''  =  ,^(A). 

Suppose  that  there  is  a  large  enough  Ai„  among  {A»}^i  ,  then  n,n  is  approxi- 
mately equal  to  ?i(A»,)  .  Thus  is  ^(Ai„) ,  a  known  quantity,  in  the  numerical 
examples  below.  When  the  data  are  noiseless,  peirameters  n„_i,rt^_2  and 
can  be  determined  by  an  asymptotic  procedure  of  the  peeling  oS  type  [3, 
Appendix  7].  However,  this  asymptotic  approach  is  not  easy  to  use  because  of  the 
rapid  accumulation  of  round-oflF  errors  and  its  sensitivity  to  noise. 

It  is  obvious  that  one  needs  some  a  priori  assumptions  about  parameters 
{ny}^!  and  {^y}^!  .  There  are  three  typical  cases  with  noise-free  synthetic 
data  or  noisy  data  for  which  the  recovery  may  be  difficult  to  carry  out.  The  first 
case  is  when  some  layers  are  very  deep,  the  second  is  when  some  layers  are  very 
thin,  and  the  third  is  when  two  adjacent  layers  have  ahnost  the  same  parameters 
Tij  .  Hence  the  numerical  examples  are  mainly  focused  on  these  difficult  cases,  yet 
a  few  intermediate  cases  are  also  presented  for  comparison. 

The  results  obtained  are  presented  in  six  tables.  Tables  1,  3  and  5  are  for  the 
two,  three  and  four  layers  with  noise-free  data.  Tables  2,  4  and  6  are  for  two,  three 
and  four  layers  with  different  noise  levels  in  the  data.  The  results  show  that  the 
algorithm  is  applicable  to  a  large  range  of  layered  structures  and  is  efficient.  In  the 
numerical  examples  the  quantities  ny  ,  Ay  ,  and  z,  are  dimensionless.  For  example 
rij  =  2  means  that  one  has  chosen  some  unit  of  and  ny  =  2  x  (unit  of  ny )  . 
For  a  practical  problem,  if  one  has  a  scale  system  under  which  the  value  of  each 
variable  falls  into  the  range  we  describe  below,  then  one  can  apply  the  method  to 
this  problem. 

First,  if  the  thicknesses  of  layers  are  not  too  small  and  the  depths  of  layers  are 
not  too  large,  then  recovery  is  reasonably  accurate.  Specifically,  we  require  that 
the  thickness  of  each  layer  should  not  be  less  than  0.2  and  the  depth  of  each  layer 
should  not  be  more  than  20 .  Under  these  assumptions,  the  recovery  is  accurate 
even  for  noisy  data,  the  majority  of  which  have  noise/signal  ratio  less  than  0.01 . 
The  case  when  layers  are  thin  and  deep  is  more  difficult  than  that  when  layers  are 
thin  and  shallow.  Also  one  will  not  have  much  difficulty  to  recover  the  parameters 
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{"jj-^^j}  if  there  is  no  big  difference  in  the  thicknesses  of  layers  even  though  they 
are  very  thin  or  they  are  deep  down  in  the  earth. 

Secondly,  although  one  can  divide  both  sides  of  (5)  by  any  number  to  make 
{^j}^!  larger,  the  level  of  noise  in  ^(A)  is  changed  at  the  same  time.  Therefore, 
the  relative  error  of  recovery  remains  basically  the  same.  Clearly,  the  larger  ny  , 
J  =  1»  2, . . . ,  m ,  the  easier  the  recovery,  provided  that  the  level  of  noise  is  the 
same.  In  fact,  the  range  of  ny  for  which  the  recovery  is  feasible  depends  highly 
on  the  difference  between  Zj  and  zy+i .  We  suggest  that  ny  be  larger  than  1 , 
J  =  1»  2, . . . ,  m  ,  if  the  level  of  noise  is  not  greater  than  0.001 .  Some  restriction  one 
should  impose  on  the  difference  of  the  parameters  of  the  adj  acent  layers:  |  n j+ 1  -  | 
should  be  of  order  100c  where  €  is  the  absolute  error  of  the  data.  No  restrictions 
are  needed  on  |n^  -  n,| ,  t  #  j  -  1,  j  ,  j  +      i,  j  =  l,2,...,m. 

Thirdly,  the  sample  rate  AA» ,  t  =  1,2, . . . ,  ^■  for  ^(A)  can  be  chosen  be- 
tween 0.01  and  0.1 .  The  number  of  samples  can  be  around  100  ,  and  this  num- 
ber should  be  increased  as  the  number  of  layers  increases.  One  should  not  choose 
large  A  as  sample  points  except  one  point  to  get  and  does  not  have  to  choose 
equidistant  samples.  The  rate  of  convergence  of  i^*)  depends  on  the  choice  of 
samples. 

Fourthly,  since  there  are  many  local  minimizers,  we  may  not  be  able  to  get  a 
global  one  if  we  start  with  an  arbitrary  initial  point.  We  recommend  to  use  the 
criterion 

{F{x')<e^    and     i|  (J(x*))^(/(x*))  ||<  e}  (lO) 

to  decide  which  minimizer  is  the  global  one.  In  order  to  get  better  initial  points 
one  can  use  an  auxiliary  optimization  method,  such  as  pattern  search  method  or 
simplex  method  [7],  and  then  apply  the  above  algorithm  to  the  prospective  points 
found  by  the  auxiliary  method. 

Fifthly,  since  data  are  always  noisy  in  practice,  it  is  important  to  show  that 
the  method  can  handle  the  cases  with  noisy  data  when  the  noise  level  is  within  a 
certain  level.  Note  that  4>{X)  decreases  very  fast  as  A  grows.  In  our  examples 
noise/signal  ratio  for  88%  of  the  data  with  noise  level  0.001  is  <  0.01  .  The 
remaining  data,  though  very  noisy,  do  not  play  a  significant  role  in  the  recovery. 
Our  results  show  that  the  algorithm  works  well  in  the  presense  of  noise.  In  tables 
2,  4  and  6,  one  can  see  that  if  the  level  of  noise  is  increased  a  little,  then  the 
accuracy  of  the  recovered  parameters  does  not  decrease  much.  At  the  recovered 
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values  of  {ny.zy}  every  fi  in  the  target  function  F  is  not  greater  than  the  noise. 
Therefore  the  accuracy  of  the  recovery  can  not  be  improved. 

We  conclude  that  the  algorithm  can  be  used  for  practical  recovery  of  the 
layered  medium  from  the  seismic  data. 

rv.  Additional  remarks.. 

The  algorithm  can  be  generalized  for  the  case  when  ny  and  Zj  are  complex 
numbers.  The  number  of  the  layers  can  be  increased  to  10  by  the  given  method 
according  to  the  results  of  numerical  experiments. 

Let  us  give  a  method  for  estimating  the  number  of  layers  from  the  data. 

Let  T(A)  :=  nje'^''  .   Apparently,  formula  (5)  and  r(A)  have  the 

same  form.  Set  =  T{k  ■  AA)  ,  e,"  =  e"^*  '^  ,  AA  >  0  ,  j  =  1, 2, . . . ,  m  ,  thus 
Tk  =  E]?Li«j^y  .  *  =  1,2,....  Assume  ^.  ^  for  t  y,  ij  =  1,2,.. .,m, 
then  we  have 

rank  Aj,  =  min(m,p)       for  ny  >  0,  j  =  1, 2, . . . ,  m,  (11) 

and 

rank  Ap  =  m       for  p  >  m.and  ny  ^  0,  j  =  1,2, . . . ,  m.  (12) 
where  rank  A  is  the  rank  of  the  matrix  A  ,  and 


/  To      Ti...     rp_i  >y 
Ti  Tj... 


VTp.i    Tp...  Tjp.j/ 


Hence  the  number  m  is  the  smallest  p  starting  from  which  rank  Ap  does  not 
change  as  p  grows.  Since  each  ^y  ,  j  =  1, 2, . . . ,  m ,  exponentially  decreases  when 
k  grows,  Ap  becomes  ill-conditioned  as  p  increases.  Therefore  it  is  difficult  to 
compute  its  rank.  Prony's  method  can  also  be  used  for  determining  ny  and  Zj  , 
J  —  1, 2, . . . ,  m  .  Prony's  method  is  compUcated  and  very  sensitive  to  the  noise  in 
the  data.  An  extensive  bibUography  on  Prony's  method  can  be  found  in  [8]  and 
[3,  Appendix  7]. 

One  can  use  (11)  or  (12)  for  estimating  m ,  and  then  use  our  algorithm  with 
the  criterion  (10)  for  recovery  of  ny  and  Zj  . 
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APPENDIX  1 


Theorem  l..  If       >  qo  :=  w'max,  |n(z)| ,  then  the  integral  equation  (2)  is 
uniquely  solvable  by  iterations,  and  u  is  analytic  in  X  in  the  region  ReX  > 
PROOF:  Let  Vii  :=  /+~  ^(A,  z  -  z')q{z')u{X,  z')dz' .  One  has 


1/2 


sup,  \g{X,z  -  z')q{z')\dz'  <  sup,  ^  f^'" 

?o  2  ?o 

By  Young's  inequality,  V  maps  i'(-oo,+oo)  into  itself,  and  <  ^  .  If 

9o  <  A'  ,  then  |1V||  <  1 .  Therefore  the  series 


u(A,  z)  =  9iX,  z)  +  Vg{X,  z)  +  ...  +  V^g{X,  z)  + 


(13) 


gives  the  unique  solution  to  equation  (2)  and  converges  uniformly  in  z  and 

1/2 

A  for  A  >  go  •  Since  for  all  positive  integers  k,  the  function  V''g{X,z)  = 
/-<^  ~ — 2A — y''~^g{\z')dz'  is  analytic  in  A  and  the  series  (13)  converges 

uniformly  in  the  region  A  >  g^'  ,  it  follows  that  u  is  analytic  in  A  if  A  >  ReX  > 
i/a 


In  the  following  theorem,  we  use  the  same  notations  as  in  the  section  IV. 

THEOREM  2 . .  Assume  that  ^  for  i  ^  j,  i,j  =  l,2,...,m,  then  (11)  and 
(12)  bold. 

We  consider  (11)  Grst.  It  is  sufficient  to  prove  that  det  Ap  =  0  if  p  >  m  and 
that  det  Ap  ^0  if  p  <  m  . 

Suppose  p  >  m,  we  want  to  show  that  det  Ap  =  0 .  Note  that  det  Ap  can 
be  represented  as  a  sum  of  terms 


"ti  lis  •  •  • »»»,  det 


(  1    ^i.---  CM 


2p-2 


(14) 


13 


where  =  1, 2, . . . ,  m,  A  =  1, 2, . . . ,  p .  Since  p  >  m ,  at  least  two  columns  of  the 
above  matrix  are  limearly  dependent.  So  the  determinant  of  the  matrix  is  zero. 
This  implies  that  det  Ap  =  0  . 

Now  we  prove  that  det  Ap  ^  0  if  p  <  m .  Consider  all  the  terms  of  the 
form  (14).  We  divide  them  into  several  groups.  In  each  group  the  p  -  tuple 
r*i>«2i  runs  through  all  the  permutations  of  the  integers  «i,tj,  ...,tp.  We 

pick  up  one  of  these  groups  and  sum  all  the  terms  in  it.  Without  loss  of  generality, 
we  assume  that  for  this  group  f»i,  ij,  ...,ip)  is  a  permutation  of  (1,  2,  p). 
Therefore  the  sum  is 


det 


ynun^^r  T^u^ie^.:  zu^^er'J 


(  1 


=  det 


1. 
6 


1  \ 


det 


p-i 


\np    rip^p...  np^-'^^J 


=  Tiifij  . . .  Hp  det 


1. 
6 


1  \' 


/  1 

ver'  e?  --  ep'-'J 


>  0. 


Since  ^»  for  i  #  j,  i,  j  =  1, 2, . . . ,  m ,  and  p  <  m  ,  we  infer  that  det  Ap  >  0 
if  p  <  m.  Therefore  we  conclude  that  rank  Ap  =  min(m,p)  under  the  assumption 
(11). 

From  the  above  argument,  it  is  clear  that  det  Ap  =  0  for  p  >  m ,  regardless 
of  the  values  of  tij  .  If  p  =  m,  then  the  above  argument  yields 
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^ 


det  Am  =  "ina  . . .  det 


(    1       1...        1  \ 

^1    f  2  •  •  •  if 


>  0. 


pTOYided  tiat  0 ,  and      7^  $y  for  i  ^  j,  i,j  =  1,2,. . .  ,m .  Therefore 

follows.  This  completes  the  proof  of  theorem  2. 
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APPENDIX  2 


In  the  section  III,  we  point  out  that  the  asymptotic  approach  practically  is 
hard  to  carry  out.  Here  we  give  an  explanation.  Consider  two  layer  case,  then  (5) 
becomes 

nj  +  (ni  -  nj)c-'^'»  -  me"'^''  =  VCA)  (15) 

Let  us  first  describe  the  asymptotic  approach  to  get  nj,zj,ni  .  Below  *  denotes 
the  recovered  parameters.  As  mentioned  before,  we  use  a  datum  ^'(•^io)  to  get  , 
where  Aj„  is  a  large  number.  Suppose  ^(A„)  is  a  datum  such  that  nie~^^'i'^  « 
(ni  -  nj)e~'*«i'=  and  I  is  small  enough.  Since  d  >     ,  we  can  have  such 

a  datum.  Then  we  approximately  have 

ln\ni  -  njl  -  2A<jZa  =  /nlV-CA^J  -  n\\ 

In  order  to  get  zj  ,  we  eliminate  the  small  quantity  I  hence 

/n|V.(A,J  -  n;| 
^2-  2A:^   (16) 

From  (14),  using  another  datum  V'(Ata) ,  we  get 

^         e-2A<3Zj  -  e-2Aj,<i 

Assume  now  that  ni  =  4,  nj  =  7,  d  =  5,  zj  =  2  .  Let  us  analyse  more  carefully 
the  way  we  get  z\  .  (16)  is  more  precise  when  A  becomes  larger.  For  example. 


,  .        /n(1.2745063  *  10-^^) 

A  =  lU,     Zj  =  —  '-  =  1 


20 


9450 


X     on       .        /n(5.4145542  *  10-"^) 

A  =  20,    z;  =  i  L  =  1.9725 

40 

One  sees  that  if  A  increases,  the  numbers  under  logarithm  decrease  extremely  feist, 
while  the  accuracy  of  z,  does  not  change  much.  Practically,  the  round-ofiF  error  of 
computers  does  not  allow  us  to  increases  A  very  much.  Moreover  we  sdways  have 
noise  in  V'(A,, )  -  nj  which  makes  the  estimation  of  zJ  difficult.  Therefore  we  are 
not  able  to  recover  zJ  by  this  approach  even  in  this  simplest  two  layer  case. 
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APPENDIX  3 


Here  is  the  list  of  arguments  of  our  program  which  is  appended. 

AA  -  Real  work  matrix  of  order  (2Ar  -  2)  x  {2k  -  2) . 
AJ  -  Real  work  matrix  of  order  (MO  +  3A;  -  3)  x  (2Jb  -  2) . 
AJF  -  Real  work  vector  of  length  {2k  -  2) . 
AJTJ  -  Real  work  matrix  of  order  (2A;  -  2)  x  (2*  -  2) . 
AJTJJl  -  Real  work  matrix  of  order  (2*  -  2)  x  (2A:  -  2) . 
ALAMDA  -  Input  real  vector  of  length  MO  ,  which  contains  the  points  where  the  data 
are  picked  up. 

AN  -  Real  work  vector  of  length  k.  In  the  end,  it  contains  the  estimation  of 
numbers  in  ANE. 

ANE  -  Input  real  vector  of  length  A: ,  in  which  the  exact  parameters  rij  ,  j  = 
1,2,...  ,k ,  are  stored. 
AZ  -  Real  work  vector  of  length  k.  In  the  end,  it  contains  the  estimation  of 
numbers  in  AZE. 

AZE  -  Input  real  vector  of  length  k  ,  in  which  the  exjict  parameters  zy  ,  j  = 
1,2,...  ,k  ,  are  stored. 
EPSN  -  Input  real  number  which  is  the  error  of  the  data. 
F  -  Real  work  vector  of  length  MO  +  3Jfe  -  3  . 
FE  -  Real  work  vector  of  length  MO  ,  which  contains  the  exact  data. 
FIl  -  Real  number  which  is  the  value  of  the  target  function. 

K  -  Input  integer  which  is  the  number  of  the  layers. 
MO  -  Input  integer  which  is  the  number  of  the  data. 
T  -  Real  work  number  which  is  the  maximum  value  among  all  the  components 
of  the  vector  {J {x^'''^))'^ f  {x^-''^)  in  formula  (7)  in  each  step  of  iteration. 
P  -  Real  work  vector  of  length  {2k  -  2)  . 
X  -  Real  work  vector  of  length  (2A;  -  2)  . 
XI  -  Real  work  vector  of  length  (2A;  —  2)  . 
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PILEI    T  FORTRAN      A        KANSAS   STATE   UNIVERSITY    VN/SP  CHS 

THE  PROGRAM  FOR  SOLVING  THE  INVERSE  LAYER  PROBLEM 
THE  MAIN  PROGRAM 

REALMS   AN(4),AI(4), ANE<*),AZE(4),AJ(89,6).AJTJ<6,6), ALAMDACaO), 
*■  X(6)fXl(6),F(89),FEC80),AJFC6),P(6),AJTJl(6,6) 

REAL*!   ALPHA.CAMNA, R,OE lEP.EPS, EPSN.FIl ♦F1,T|AA(6,12),PHI 
DATA   ALAHOA/0.02O*O0t  0. 040*00 i0.060*00«0.080*00>0.100«00i 
«  0.120*00,0. 140*00, 0.  160«00t  0.  180*00.0. 200*00  . 

•  0.21 0*00, 0.220*00, 0.230*00, 0.240*00. 0.250*00. 

*  0. 260*00, 0.2TO*00, 0.280*00, 0.290*00, 0.300*00, 
«  0.310*00, 0.320*00, 0.330*00. 0.340*00, 0.350*00. 

•  0.360*00, 0.370*00, 0.390*00, 0.39  0* 00, 0.400*00, 

*  0.420*00,0.4*0*00,0.460*00,0.480*00,0.400*00, 

•  0.520*00, 0.540*00, 0.560*00. 0.580*00. 0.500*00. 

•  0.620*00,0.640*00,0.660*00,0.680*00.0.700*00. 

•  0.720*00, 0.740*00, 0.760*00, 0.730*00, 0.800*00, 

*  0.820*00, 0.940*00, 0.360*00, 0.830*00, 0.900*00, 
«  0.920*00,0. 940*00, 0.960*00, 0.980*00. 1.000*00, 

*  1.020*00, 1.040*00, 1.060*00. 1.080*00.1. 000*00, 
«  1. 120*00, 1. 140*00, 1.160*00,1. 180*00, 1.200*00, 

•  I. 220*00, 1.24 0*00, 1.260*00, 1.280*00. 1.300*00, 

•  1. 320*00, 1. 340*00, 1.360*00,1. 380*00, 1.400*00/ 
DATA  ANE/5. 00*00, 8. 00*00, 3. 00*00, 7. 00*00/ 

DATA  A2E/6. 00*00, 4. 50*00, 3. 10*00, 1.50*0  0/ 

DATA  X/4. 98 100*00,  7. 8 20*00, 2.  ?.60*00, 4. 540*00, 3.0 70* 00, 1.5010*03/ 
K«4 

M32«K-2 

NN«2*N 

M0"80 

M»M0*3»K-3 
ALPHA*1 . 00*00 
GAMMA15 . 00*00 
EP»1. 00-09 
EPS'l. 00-09 
EPSN»0. 00010*00 
R=500. 00*00 
00   2   1=1, H-1 

ANcn.xci) 

2  AZ(I*1)=XCI*H-1) 

AN(K)eANE(K)*0. 10-09 

Ai(i)«4ze(i) 

00   10    1  =  1, 10 

CALL    FF( AHE , AZE .H, ALAMOACIi ,PHI) 
WRITE(6,8)    I, PHI 

8  FOi»MAT(2X,2HI  =  ,I6,4X,3MFE-,02  0.  10) 

10         FE( I ) »PHI*EPSN 

WRIT£(6,20)  ANE,AZE 
2  0  FORMATC1X.3HN1.,014.6.2X. 3HN2.,C14.6, 

♦  2X.3HN3»,014.6,2X,3MN4=,D14.4/ 

♦  1K,3HZ1»,D14.6,2X,3HZ2=.D14.6, 

♦  2X.3HZ3«.014.6,2X,3HZ4=,014.6) 
HRITE(6,24)    AN, 42 


It 


'  ^     •  A       KANSAS   STATE  UNIVERSITT   VN/SP  CHS 


2*         FORHATdK.lSHINITIAL  VALUE/ 1 X  ,  3MN 1  -  ,  01     6  ,  3X  ,  3HN2- ,  01*.  6  ,  3X  , 
*  3MN3-t014.6, 3X,3HN4«,014.6/ 

00  40  I-l.NO 

CALL  FFCAN, A2,K,ALAN0A(I},PHI) 

40  F{I)»PHI-F£<I) 

00   42  I>H0*1.N 
♦2         F(I)-0. 00*00 

00   50   1*1, N 

IF(X(I))    44, 50, 50 

F(MO*I)»R«X(I)»X(I) 
SO  CONTINUE 

00   56  I«1,X-1 

IF(4Z(I)-a:(I»i)>  S4.56.S6 

5*  ''<«0*N*I).(AI(I)-AZ(I*n)«<AZ(I)-A2Cl»n)*« 

56  CONTINUE 

60         WRIT6(6,iO)  AN,A2 

FI-0. 00*00 

00   70  I<1,M 
70  FI»FI*F(I)«FCI) 

FI«OSQRT<FI/OFlO*T(H)) 

MRITEC6,72)  Fl 

72  FORMATC IX,5MFSUM-,016.8) 

30  IW<-1 

IF( ALPHA. LT. 0.0000050*00)  GOTO  79 
ALPHA»ALPHA/GAHHA 

79  CALL  «J(AN,Ai,ll,AJ,N,l«,MO,ALAHOA,R) 
00  100   I«1,N  ^ 

00  80  I1-1,N 

80  AJTJ(I,n)>0. 00*00 
AJF(I)«0. 00*00 

100       P(I)«0. 00*00 

DO  UO  I>l,N 

00   no  I1>1,N 

DO  110  I2»1,H 
XIO  *JTJ(I,Il)»*JTJ(I,n)*Aj<i2,i)«4j(i2 

00   120  I3»1,M 
120  *JFCI).AJF(I)*4j<l3,i),F(i3, 
140  CONTINUE 

00   160  I-1,N 

DO   160  ll«l,M 
160  AJ7J1(I,I1)-»JTJ<I,H) 
170       WRITE<6,142)  AJF 
1*2  F0RMATC1X,4HAJF=,3020.10) 

00   180  I<1,N 

00   180  J>1,N 
180       AJTJCI, J).AJTJ1(I,J> 

00   182  I>1,N 
182  *JTJ(I,I).AJTJ(I,I)*ALPHA 

IFCALPHA.LT. 0.010*00)   GOTO  185 

00  184  I'l.N 

DO   184  Jsl,N 

184  *JTJ(I, J).AJTJ<I,J)/ALPHA 

185  CALL  IVSN(AJTJ,AA,N,NN,0E,EP,IS) 
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'^^^'ri*'-  T  FORTRAN     «       KANSAS  STATE  UNIVERSITT  VM/$P  CMS 

IT"IT4l 

IF(IS.CT.O)  GOTO  200 
WRITE(6tl90)  IS 

190       F0RM*I(1X,25HTME  HATRIX   IS  OECENE R ATEO, I 1 0 > 

GOTO  230 
200       00  202  I'l.N 
202       P(I)«0. 00*00 
00  210  I>1,N 
00  210  I1>1.N 

210  PCO»P(I)»AJTJ(r,Il)«AJF(Il) 
hRITE(6,10S}  P 

105  FORMAT(IX,2HPP,3020.10) 
URITE(6,106)  ALPHA 

106  F0fiHAT(lX,6HALPHA«,015.8) 
WRITE(6i84)  IT 

84  F0RMAT(5X,3HIT-,I6) 

211  I>1 

212  IF(P<I)-1. 00*00)   213, 214, 214 

213  I>I*1 
IF(I.GT.N)  GOTO  219 
GOTO  212 

214  00  21S   1*1, N 

215  PC  I  )>!P(  I  5/2.00*00 
GOTO  211 

219  00  220  I>1,N 

220  X1<I)«X(I)-P(I} 
00   222    I-l  ,»  -1 

ANcn-xun 
222  Ai(i»n«xui.n-n 

00  230  I>1,H0 

CALL  FF(AN,AZ,H,ALAMOA<I),PMI) 

230  F(I)»PMI-Fg(I) 
00   2J1  I>M0*1,M 

231  F(I)sO. 00*00 
00  240  I'l.N 
IF(X1(I))  238,240,240 

238  F(M0*I).XUn*xi(n«R 
2  40  CONTINUE 

00   246  I>1,K-1 

IF(«an-42(lM)>  244,246,246 

24^  ^j;;;j;j^'^»^«i>-«i(i»»>>*<*Ki)-Az<i*i),»R 

249  MRITE(6,20)  AN.A2 
FI1«0.0 
00  2S2  I>1,H 

252  FIl.Fn«F(I)«F<I) 

FI1=0SQRT(FI1/0FL0ATCM)) 
MRITE{6,72>  fU 
IFCFIl.LT.FO   GOTO  270 

250  lu<lw*l 

ALPHA>ALPHA*GAMHA 
HRITE(6,260)  lU 
260  F0RHAT(1X,3MIM.,15) 
T«0 

00   264  I-l.N 
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FORTRAN 


I  • 


264 


266 


270 
280 


284 


290 


999 


10 


12 


20 

26 
28 
30 


«0 


50 


130 


K*MS»S  STATE   UNIVERSITT   VM/SP  CHS 


IF(0A8$(«JF(I)).LE.T)   GOTO  264 

T"0ABS(AJF(1)) 

CONTINUE 

MRITE(6.266>  T 

FORMAT(1X,5HHAJF»,020.10) 

IFCT.LT.EPS)  GOTO  999 

IFCIM.GT.T)  GOTO  999 

GOTO  170 

URITE(6.280)'  ALPHA 
F0RHAT(1X,6HALPHA>,014.8) 
IF(IT.GT.IOO)   GOTO  999 
T«0 

00  284  I«1,N 

1F(DABS<AJF(I)).LE.T)  GOTO  284 

T-0A8SCAJF(I)) 

CONTINUE 

MRITE(6,266)  T 

IFCT.LT.EPS)  GOTO  999 

00  290  I>1,N 

XCn-XH  I) 

FI-Fli 

GOTO  30 

STOP 

ENO 

R£AL»8  ANCK),A2<K),AJ<M,N),ALAHDA(M0> 

REAL«8  R 

00   10  J>1,K-1 

DO   10  I>1,N0 

AJ(I , J).0£XP(-2.00»00«ALAMOA( I)»AIC  J«l)) 
-OEXP(-2.00*00*ALAMOA(I)»»i( J)) 

AJ(  I,  J»K-n«2.00»00»{AN(  J»l  )-AN(  J))«ALAMDA(n 

»OEXP(-2.00*00*ALAMOA(I)«AZ( J*l)) 
00   12  I<M0*1,M 
00   12  J«1.N 
AJ(I,J}>0. 00*00 
00    30  I«l,K-l 
IFCAN(I))  20,26,24 
AJ(H0*I,I)»2.0O*0O»AN<I)«R 
IFCAZCI*!))  28,30,30 

AJ(H0*K-1  ♦I,K-l*I).2.0O*00«A2(I*n«R 

CONTINUE 

00   50  I»1,H-1 

IF(AI(I)-AZ(IM))   40,  50,  50 

AJ<«0»NM.*-2»I).-2.00»00«(AZ(I*l)-*l(l))«R 

AJCMO*N»I ,n-l*I ).-4j(M*N»I,K-2*I) 

CONTINUE 

AJ(M0»N»l,H-l)«O. 00*00 
RETURN 
ENO 


SUdROUTINE  FF(ANS,AIS,K,C,PHI) 
REAL«8  ANS(K),A2 
R£AL*8  C.PMI 


* 

Zi 


file:  T 


FORTRAN 


K*NS*S   STATE  UNIVERSITT   VM/SP  CMS 


PHI-0. 00*00 

00   10  I-ltK-1 

10         PHI»PMI**MS<I)*<0EXP(-1.00»00*C«AZS(1*1  )) 

•  -0EXP(-2.00«00*C«AZS(I})) 

PHI-PMI*AMS<K>»<1.00*00-OEXP{-2.00*00«C*AZS<K))) 
RETURN 

ENO 

SUBROUTINE  IVSN(A,AA,N,NNtOE.EP,IS) 
REAL*8   A(N>N) lAACN.NN) 
REAL«8  T.OE.EP.C 
IS-1 

OE>1.00«00 

00  6  I-l.N 

00   7  J«X,N 
7  AACI,  J)»A(I,J) 

00  8  J-liN 
e  AA(I, J*N)«0.0 

6  AA(  I,  I*N}-1 .0 

00   100  I>1,N 

T'O. 00*00 

K»I 

DO    20    J>I ,N 

IF(0AB$CAACJ,I)).L£.OABSCy)«1.0O-10)  GOTO  20 

T»AA( J, I > 

20  CONTINUE 

OE>oe»r 

IF(OABS(T).CT.EP)  GOTO  30 
IS«-1 
RETURN 

30      iF(K-n  4o,soi4o 

40  00   4*  J»I,NN 

C=AA(I, J) 

AA(  I, J)aAA(K, J)/T 
♦«  AA(R,J).C 

DE»-OE 

GOTO  62 
SO         00  60  J>I,NN 
60  AA( I , J)=AA( I, J)/T 

62         00  80  J23l*l,2*N 

00  80  J1>I*1,N 
80         *»C Jl. J2).AA(J1,J2)-A»<J1,I)«AA(I,J2) 
100  CONTINUE 

00   120  I-1,N-1 

00   120  J2:N*1,NN 

00   120  J1>1,N-I 

120  *«<JI.J2J«AA(J1,J2)-AA<J1,N»1-1)»AA(N*1-I,J2) 

00    140  I.l.N 

00   140  j3l,N 
1*0        A(I , J).AA(I, J»N) 

RETURN 

ENO 


TABLE  1:  TWO  LAYERS,  NOISELESS  DATA  ♦ 


"1 

4.0000 

12.0000 

7.2000 

12.0000 

4.0000 

15  0000 

15  noon 

112 

2.0000 

7.0000 

7.0000 

7.0000 

7.0000 

7.0000 

zi  (given) 

5.0000 

30.0000 

14.0000 

14.0000 

19.0000 

5  0000 

3.0000 

24.0000 

8.0000 

13.0000 

18.0000 

4.2000 

4.8000 

4.0001 

11.9998 

7.2000 

11.9800 

4.1029 

15.0000 

15.0742 

2.0000 

7.0000 

7.0000 

7.0000 

7.0000 

7.0000 

7.0000 

^2 

3.0000 

23.9993 

8.0000 

12.9962 

17.9654 

4.2000 

4.8018 

t  all  the  results  are  rounded  to  four  decimals. 
*      denotes  the  recovered  parameters. 
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TABLE  2:  TWO  LAYERS,  NOISY  DATA 


"1 

4.0000 

12.0000 

7.2000 

15  0000 

2.0000 

7.0000 

7.0000 

7.0000 

zi  (given) 

5.0000 

30.0000 

14.0000 

5.0000 

3.0000 

24.0000 

8.0000 

4.6000 

"I 

3.9990 
3.9896 

11.9953 
11.9543 

7.2003 

15.0351 
15.3714 

noise  level  0.0001  ^ 
noise  level  0.001 

"1 

2.0001 
2.0010 

7.0001 
7.0010 

7.0001 

7.0001 
7.0010 

noise  level  0.0001 
noise  level  0.001 

^5 

2.9987 
2.9871 

23.9934 
23.9354 

7.9929 

4.6016 
4.6146 

noise  level  0.0001 
noise  level  0.001 

t     by  noise  level  we  mean  the  absolute  error  of  the  data  in  tables  2,  4  and  6. 
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TABLE  3:  THREE  LAYERS,  NOISELESS  DATA 


10.0000 

10.0000 

4.0000 

4.0000 

8.0000 

8.0000 

7.0000 

7.0000 

5.0000 

5.0000 

2.0000 

2.0000 

4.0000 

4.0000 

1.0000 

1.0000 

5.0000 

5.0000 

zi  (given) 

3.0000 

20.0000 

4.0000 

4.0000 

0.9000 

4.0000 

2.0000 

15.0000 

3.8000 

1.2000 

0.6000 

3.7000 

1.0000 

13.0000 

1.0000 

1.0000 

0.3000 

3.4000 

«I 

10.0000 

9.9978 

4.0791 

4.0000 

7.9998 

8.0505 

7.0000 

6.9404 

5.0003 

4.9649 

1.9997 

3.0246 

4.0000 

4.0000 

1.0000 

1.0000 

5.0000 

5.0000 

2.0000 

14.9726 

3.7823 

1.2043 

0.5998 

3.7375 

1.0000 

12.9841 

1.0000 

0.9992 

0.3001 

3.3320 

2S 


TABLE  4:  THREE  LAYERS,  NOISY  DATA 


"1 

10.0000 

4.0000 

4.0000 

8.0000 

7.0000 

5.0000 

5.0000 

2.0000 

4.0000 

1.0000 

1.0000 

5.0000 

zi  (given) 

20.0000 

4.0000 

4.0000 

0.9000 



15.0000 

3.8000 

1.2000 

0.6000 

"3 

13.0000 

1.0000 

1.0000 

0.3000 

"J 

10.0307 
10.4230 

4.2136 
4.5122 

4.0001 
4.0053 

8.0263 

noise  level  0  0001 
noise  level  0.001 

7.2731 
7.6933 

5.0002 
4.9955 

4.7976 
4.7888 

2.0486 

noise  level  0  Dnm 
noise  level  0.001 

4.0001 
4.0010 

1.0001 
1.0010 

1.0001 
1.0010 

5.0001 

noise  level  0.0001 
noise  level  0.001 

^2 

15.1933 
16.2958 

3.7468 
3.6303 

1.2261 
1.1991 

0.6029 

noise  level  0.0001 
noise  level  0.001 

4 

13.0488 
12.9653 

1.0000 
0.9984 

0.9947 
0.9912 

0.2985 

noise  level  0.0001 
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TABLE  5:  FOUR  LAYERS,  NOISELESS  DATA 


5.0000 

5.0000 

5.0000 

5.0000 

5.0000 

5.0000 

8.0000 

8.0000 

8.0000 

8.0000 

8.0000 

8.0000 

ns 

6.0000 

6.0000 

6.0000 

6.0000 

6.0000 

6.0000 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

zi  (given) 

5.2000 

7.0000 

7.5000 

7.5000 

8.0000 

4.5000 

3.8000 

6.2000 

2.4000 

6.4000 

7.1000 

3.3000 

2.5000 

1.8000 

1.7000 

5.3000 

6.2000 

2.3000 

z* 

1.2000 

1.0000 

1.0000 

1.0000 

5.2000 

1.1000 

nl 

5.0293 

5.0018 

5.0005 

5.0062 

5.0106 

4.9584 

8.1760 

8.0001 

8.3189 

8.0256 

7.4764 

7.7602 

^3 

6.0146 

6.0004 

6.0622 

6.0000 

5.6168 

5.9792 

< 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

3.7460 

6.1994 

2.3692 

6.3920 

7.2012 

3.3662 

2.5531 

1.8002 

1.7692 

5.3054 

5.6579 

2.2304 

< 

1.1991 

1.0000 

0.9973 

1.0000 

5.1973 

1.1012 
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TABLE  6:  FOUR  LAYERS,  NOISY  DATA 


5.0000 

5.0000 

5.0000 

5.0000 

5.0000 

5.0000 

"2 

8.0000 

8.0000 

8.0000 

8.0000 

8.0000 

8.0000 

"3 

6.0000 

6.0000 

6.0000 

6.0000 

6.0000 

6.0000 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

zi  (given) 

5.2000 

7.0000 

7.5000 

7.5000 

8.0000 

4.5000 

3.8000 

6.2000 

2.4000 

6.4000 

7.1000 

3.3000 

«3 

2.5000 

1.8000 

1.7000 

5.3000 

6.0000 

2.3000 

1.2000 

1.0000 

1.0000 

1.0000 

5.2000 

1.1000 

—  * 

5.0369 

5.0430 

5.0011 

5.0021 

5.2246 

5.0721 

noise  level  0.0001 

5.1971 

6.1799 

noise  level  0.001,  0.005+ 

'"2 

8.1152 

8.0007 

7  77S9 

tilt  0£i 

7  171R 

r  .0400 

noise  level  0.0001 

8.3570 

7.4567 

noise  level  0.001,  0.005 

"a 

5.9956 

5.9911 

6.0363 

5.9997 

5.6422 

5.9527 

noise  level  0.0001 

5.8726 

5.9922 

noise  level  0.001,  0.005 

< 

10.0001 

10.0001 

10.0001 

10.0001 

10.0001 

10.0001 

noise  level  0.0001 

0.0010 

10.0050 

noise  level  0.001,  0.005 

3.7520 

6.1881 

2.3657 

6.4472 

7.1332 

3.2890 

noise  level  0.0001 

3.5874 

5.0568 

noise  level  0.001,  0.005 

2.5229 

1.7979 

1.7603 

5.2356 

5.5113 

2.2177 

noise  level  0.0001 

2.4738 

4.4969 

noise  level  0.001,  0.005 

1.2007 

1.0007 

0.9989 

1.0001 

5.1684 

1.1031 

noise  level  0.0001 

1.2127 

1.0015 

noise  level  0.001,  0.005 

+  The  data  in  the  second  column  have  noise  level  0.0001  and  0.001,  and  the  data  in  the  fifth  column 
have  noise  level  0.0001  and  0.005. 
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In  this  paper  an  inverse  problem  is  formulated  and  its  numerical  solution  is 
given.  The  problem  is  of  interest,  for  example,  in  geophysics.  Numerical  results 
show  that  the  algorithm  is  efficient  for  noisy  data.  Prjictical  recommendations  are 
given.  The  computer  code  which  solves  the  inverse  problem  is  appended. 


